Dynamical Determination of the Innermost Stable Circular Orbit of Binary Neutron 

Stars 



Pedro MarronettiQ Matthew D. Duez, and Stuart L. ShapircQ 
Department of Physics, University of Illinois at Urbana- Champaign, Urbana, IL 61801 



o 
o 

(N 

X> 

<D 

O 

(N 

(N 
> 

m 
o 

(N 

m ; 
o ■ 

cr 

i ' 



x 



Thomas W. Baumgartqj] 

Department of Physics and Astronomy, Bowdoin College, Brunswick, ME 04011 

We determine the innermost stable circular orbit (ISCO) of binary neutron stars (BNSs) by 
performing dynamical simulations in full general relativity. Evolving quasiequilibrium (QE) binaries 
that begin at different separations, we bracket the location of the ISCO by distinguishing stable 
circular orbits from unstable plunges. We study T = 2 polytropes of varying compactions in both 
corotational and irrotational equal-mass binaries. For corotational binaries we find an ISCO orbital 
angular frequency somewhat smaller than that determined by applying turning-point methods to 
QE initial data. For the irrotational binaries the initial data sequences terminate before reaching a 
turning point, but we find that the ISCO frequency is reached prior to the termination point. Our 
findings suggest that the ISCO frequency varies with compaction but does not depend strongly on 
the stellar spin. Since the observed gravitational wave signal undergoes a transition from a nearly 
periodic "chirp" to a burst at roughly twice the ISCO frequency, the measurement of its value by 
laser interferometers (e.g LIGO) will be important for determining some of the physical properties 
of the underlying stars. 

PACS numbers: 04.30.Db, 04.25. Dm, 97.80.Fk 



The emission of gravitational radiation drives the slow 
inspiral of neutron star and black hole binaries towards 
their coalescence and merger. Fully general-relativistic 
numerical simulations are required for the accurate de- 
scription of the late inspiral and plunge epochs of the bi- 
nary evolution (see, e.g., Q] for a recent review). While 
complete orbits of binary black holes have not been nu- 
merically simulated yet, simulations of binary neutron 
star (BNS) mergers are now becoming sufficiently ma- 
ture to provide results of astrophysical interest (e.g. @])- 

One piece of information of great astrophysical inter- 
est is the location of the innermost stable circular or- 
bit (ISCO). The evolution of a binary system occurs in 
three distinct phases (1) A slow, adiabatic inspiral 
phase that is driven by gravitational radiation reaction 
forces and can be approximated as a sequence of quasi- 
circular orbits; (2) a brief transition phase, where the 
inward radial motion increases and the orbital motion 
changes from slow inspiral to rapid plunge; (3) a plunge 
phase, terminating in the merger of the objects. The 
ISCO resides within the transition region; its identifica- 
tion is complicated by the fact that it is not arbitrarily 
sharp and cannot be localized precisely. The gravita- 
tional wave quasi-periodic "chirp" signal of the slow bi- 
nary inspiral ends at about the twice the orbital angular 
frequency of the ISCO, where it changes its form to a 
wave signal characteristic of a burst (compare |J|). 

Within the framework of Newtonian and post- 
Newtonian gravity, the ISCO has been determined by 
different methods (see, e.g., the reviews [HEt and refer- 
ences therein). Much less is known for fully relativistic 
binaries. For corotational binaries, a turning point on a 



curve of the binding energy vs. separation for quasiequi- 
librium (QE) models along a sequence of constant rest 
mass marks the onset of secular instability 0> • I n the 
following we will refer to this point as the QE-ISCO. No 
such theorem exists for irrotational binaries or for the 
onset of dynamical instability 0- Locating the ISCO 
dynamical instability therefore requires dynamical evo- 
lution simulations of the full set of Einstein's equations 
for the gravitational field, coupled to relativistic hydro- 
dynamics in the case of BNSs. 

In this paper we present the first attempt to dynam- 
ically locate this ISCO. We identify BNS configurations 
that correspond to stable and unstable circular orbits 
by evolving binary initial data sets for different separa- 
tions. The objective is to bracket the location of the 
ISCO by distinguishing configurations that can maintain 
quasi-circular motion for more than one orbital period 
and systems that plunge and coalesce in a fraction of 
that time. 

We adopt the QE initial data presented by Marronetti 
and Shapiro [lfl ]. describing two identical neutron stars 
in quasi-circular orbit. These data have been constructed 
using the conformal thin-sandwich (CTS) decomposition 
of the constraint equations, together with maximal slic- 
ing and spatial conformal flatness. The formalism intro- 
duced in ^3 allows for a free specification of the spin 
of the stars in an approximate fashion. In this paper we 
consider corotational binaries as well as "irrotational" bi- 
naries \Ti\ with zero (equatorial) fluid circulation, which 
are believed to be more realistic astrophysically [l2j . Se- 
quences of corotating binaries feature a minimum in the 
binding energy at the QE-ISCO. Irrotational binaries, 
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however, typically terminate before reaching this mini- 
mum (see Figs. 7 and 8 in Efl; also [HQ as well as Q 
for more references). 

We adopt a T = 2 polytropic equation of state (EOS) 
for which the maximum rest mass (gravitational mass) 
of a star in isolation in nondimcnsional units |15| is 
too = 0.180 (m = 0.164) with a compaction ratio of 
(m/R)oo = 0.216 . We study models that have two dif- 
ferent compaction ratios in isolation: a moderate value 
(to / R)oo — 0.142 (both corotational and irrotational bi- 
naries; cases A and B) and a high value (m/R) 00 = 0.195 
(only irrotational binaries; case C). These compaction 
ratios correspond to individual stars with rest masses 
too = 0.1469 and to = 0.1767, respectively. The fully 
general relativistic hydrodynamical code employed for 
this study has been introduced in Duez et al. |ljj. We 
evolve the gravitational fields using the BSSN formal- 
ism with a Courant factor of 0.46. We approximate 
maximal slicing with a "K-driver" and use a "Gamma- 
driver" shift condition that keeps f 1 = ^ lm T \ m approx- 
imately constant. The simulations were performed on 
uniform Cartesian grids in a reference frame that rotates 
with the binary, which improves conservation of angu- 
lar momentum [l|| ll3 and reduces the spurious eccen- 
tricities of the stable orbits. In all cases presented here 
the spatial volume covered by the grids is completely en- 
closed by the light cylinder (defined by coordinate radius 
Rl = 1/riorb)- A more detailed description of the evolu- 
tion of the gravitational and hydrodynamical fields, the 
boundary conditions and their numerical implementation 
can be found in Duez et al. 

Results Figures ^ 121 an d ED show the evolution of 
the coordinate separation d between maximum baryonic 
density points in each star for cases A, B and C . The 
filled circles mark the points of surface contact for those 
runs that result in a merger. In each plot we present re- 
sults for different grid resolutions and bounding box sizes, 
which are listed in Table[Q We use the highest quality re- 
sults to bracket the ISCO; these are labeled as Stable and 
Merger in the plots. Results obtained on smaller compu- 
tational grids agree with these brackets fairly well. Note 
that all three merger cases experience surface contact 
(and the related mass interchange) well after the start 
of the inspiral plunge. The ISCO parameters for each of 
the three cases are estimated as the average of the pa- 
rameters corresponding to the bracketing orbits labeled 
Stable and Merger on each figure, while the "errors" span 
the difference. We note that these "errors" are partly due 
to numerical errors (see below), and partly by the concep- 
tional difficulty of defining a sharp "ISCO" . The results 
are included in Table ITU For the corotating sequence A 
we also compare with the QE-ISCO at the QE turning- 
point. The irrotational sequences considered in this pa- 
per terminate before reaching a turning point (compare 
|l3l H^jl. The termination of a sequence indicates that 
equilibrium models do not exist at smaller separations, 
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FIG. 1: Coordinate separation vs. time for sequences A. 
The separation d is the coordinate distance between points 
of maximum rest mass density and is given in units of the 
total rest mass of the binary Mo. The curves show runs with 
different grid sizes and resolutions which are detailed in Table 
Q The filled circles mark the time of surface contact for the 
merger orbits. The empty square on the y axis marks the QE 
estimation of the ISCO separation. 




FIG. 2: Coordinate separation vs. time for sequences B. 



but since the numerical determination of this point relies 
on the breakdown of a numerical (equilibrium) code, its 
accuracy is somewhat uncertain. For both sequences we 
find that the ISCO is reached before the sequence ter- 
minates. For the irrotational sequence C one can also 
compare with the first-order post-Newtonian ellipsoidal 
results of Lombardi et al. |20j, who find an angular ve- 
locity of flmo — 0.0226 for (m/R)^ = 0.2 binaries. 

We consistently find that orbits become unstable be- 
fore the onset of instability as determined by QE meth- 
ods, meaning at a smaller orbital frequency. This is not 
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FIG. 3: Coordinate separation vs. time for sequences C. 
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TABLE I: Grid sizes and resolutions. The resolution (Res.) is 
given in number of grid points across the stellar diameter. The 
bounding box length B gives the extent of the physical space 
covered in each direction in units of total rest mass Mo (i.e.; 
the numerical grid spans from[ — B, 0.0, 0.0] to [ B, B, B] 
since we make use of the equatorial and tt symmetries of the 
systems.) 



surprising, since the orbital decay becomes fairly rapid 
just outside the QE-ISCO (compare, e.g., Fig. 2 in [j]), 
so that our criterion for merger orbits applies to those 
binaries. This result is also consistent with earlier sug- 
gestions that the transition through the ISCO may be 
fairly gradual The similarity between the corotating 
and irrotational values for the (m/R) ^ = 0.142 orbits 
suggests that the dependence of the ISCO parameters on 
the stellar spins is not strong. 

During all simulations we monitored the Hamiltonian 
and Momentum constraints as well as the conservation of 
the total ADM mass M and angular momentum J [2l| 
(the rest mass Mq is conserved identically in our evolu- 
tion scheme). An example of these for case B is shown in 
Fig. 0] In all our simulations all quantities are conserved 
very well up to merger, after which hydrodynamical ef- 
fects including shocks and shear are handled only crudely 
by our artificial viscosity scheme. Stable runs ultimately 
break down due to accumulation of numerical error. We 
find that the latter is sometimes dominated by hydrody- 



namical effects, leading to deviations in the angular mo- 
mentum, and sometimes by gravitational effects, leading 
to violations of the constraint equations. These effects 
are improved by increasing the grid resolution and the 
separation to the outer boundaries, as well as using a co- 
ordinate system that rotates with the binary as closely 
as possible. 

Summary We present prototype simulations to de- 
termine dynamically the ISCO of BNSs. Evolving QE 
initial data at different separations, we bracket the ISCO 
by distinguishing stable orbits, that remain in approxi- 
mately circular orbit for well over a period, from unstable 
ones, which decay within a period. The uncertainty in 
our results is caused both by numerical error and the 
conceptual difficulty in denning a sharp ISCO. Consis- 
tent with earlier results 0, 0] we find that binary orbits 
start to plunge somewhat outside of the "QE-ISCO" as 
determined by turning-point methods applied to QE ini- 
tial data, resulting in a cut-off in gravitational "chirp" 
signals at somewhat smaller frequency. Our preliminary 
results also seem to indicate that the dependence of the 
ISCO parameters on the stellar spins is not very strong. 

One source of error is our assumption of zero radial 
velocity in our binary initial data. More realistic ini- 
tial data at finite binary separation would incorporate 
a radial velocity that corresponds to the non-zero rate 
of inspiral at that separation. Miller 23] indicates that 
this approximation may lead to non-negligible error, es- 
pecially for black hole binaries. However, for the neutron 
star binaries and separations we consider here, the ra- 
dial velocities would be at most 1 — 3% of the tangential 
velocity 0, . Consequently, the error introduced by 
neglecting this component is likely to be smaller than 
the error bars already provided in Table ITTl 

It was recently pointed out that in dynamical evolu- 
tions of CTS initial data describing corotating BNSs, the 
four- velocity u a quickly deviates from being proportional 
to an exact helicoidal Killin g ve ctor, as assumed in con- 
structing the initial models |24J. We confirm this result 
but find that this deviation arises from a small readjust- 
ment of the gravitational fields |25|] , and not from appre- 
ciable changes in the density and velocity profiles. We 
find no evidence of a significant breakdown of quasiequi- 
librium for stable orbits and any spurious orbital eccen- 
tricities sharply decrease with increasing grid size. 
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the National Center for Supercomputing Applications at 
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port through the Fortner Fellowship at the University of 
Illinois at Urbana-Champaign (UIUC). 
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Case rotation mo (m/R)oo Mi Ji/Mf d/Mo Q mo mo (QE) }gw mi. 4 (kHz) 

~A corotating 0.1469 0142 0.2708 ± 0.0001 1.08 ± 0.01 9.0 ± 0.8 0.0162 ± 0.0021 0.0179 0697 

B irrotational 0.1469 0.142 0.2702 ± 0.0003 0.97 ± 0.04 9.4 ± 0.9 0.0154 ± 0.0022 0.662 

C irrotational 0.1767 0.195 0.3171 ± 0.0009 0.93 ± 0.04 7.7 ± 0.8 0.0199 ± 0.0029 0.838 



TABLE II: Summary of our binary cases A, B, and C and results for their ISCO. For each sequence, we show the rotation state, 
the rest mass of the individual stars mo, the compaction in isolation (m/R)^, the total initial ADM mass Mi and angular 
momentum Ji/Mf, as well as the binary coordinate separation at the ISCO d/Mo (where Mo is the total rest mass 2mo) and 
the corresponding orbital angular velocity Q mo- For the corotational sequence we also compare with the QE result for Q mo 
[To||. The wave frequency at the ISCO is few and mi. 4 is the stellar gravitational mass in units of 1.4 Mq. 
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FIG. 4: Quality control for the (rn/R)^ — 0.142 irrotational 
runs. We show from top to bottom the total gravitational 
mass, angular momentum, the Hamiltonian constraint, and 
the avera ge o f the components of the momentum constraint 
vs. time |22|. The curves correspond to the runs from Fig. 
Q labeled as Stable (solid) and Merger (dashed) . The filled 
circle marks the time of surface contact for the merger orbit. 



Fortner Fellow 

Department of Astronomy & NCSA, University of Illinois 
at Urbana-Champaign, Urbana, IL 61801 
Department of Physics, University of Illinois at Urbana- 
Champaign, Urbana, IL 61801 

T. W. Baumgarte and S. L. Shapiro, Phys. Rept. 376, 
41 (2003). 

M. Shibata and K. Uryu, Phys. Rev. D 61, 064001 (2000). 
M. Shibata and K. Uryu, Prog. Theor. Phys. 107, 265 
(2002). 

A. Ori and K. S. Thorne, Phys. Rev. D 62, 124022 (2000); 
A. Buonanno and T. Damour, Phys. Rev. D 62, 064015 

(2000) . 

M. D. Duez, T. W. Baumgarte, S. L. Shapiro, M. Shibata, 
and K. Uryu, Phys. Rev. D 65, 024016 (2001). 
T. W. Baumgarte, in Astrophysical Sources of Gravita- 
tional Radiation', edited by J.M. Centrella (AIP Press), 

(2001) . 

F. A. Rasio and S. L. Shapiro, Class. Quant. Grav. 16, 



[10] 
[11] 

[12] 
[13] 
[14] 
[15] 
[16] 
[17] 

[18] 
[19] 

[20] 
[21] 



[22] 



[23] 
[24] 



Rl (1999). 

J. L. Friedman, J. R. Ipser, and R. D. Sorkin, Astrophys. 
J. 325, 722 (1988). 

T. W. Baumgarte, G. B. Cook, M. A. Scheel, S. L. 
Shapiro and S. A. Teukolsky, Phys. Rev. D 57, 6181 
(1998). 

Approximating the stars as ellipsoids, it has been shown 
that the dynamical orbital instability occurs quite close 
to the secular instability for irrotational binaries in both 
Newtonian and first-order post-Newtonian theory (D. 
Lai, F. A. Rasio and S. L. Shapiro, Astrophys. J. Suppl. 
88, 205 (1993)). 

P. Marronetti and S. L. Shapiro, Phys. Rev. D 68, 104024 
(2003). 

For simplicity, in the reminder of the paper we will refer 
to the zero circulation models as irrotational. See [Toj| 
for comparisons between these models and irrotational 
binaries. 

C. S. Kochanek, Astrophys. J. 398, 234 (1992). L. Bild- 
sten and C. Cutler, Astrophys. J. 400, 175 (1992). 

K. Uryu, M. Shibata and Y. Eriguchi, Phys. Rev D 62, 
104015 (2000). 

E. Gourgoulhon, P. Grandclement, K. Taniguchi, J.-A. 
Marck and S. Bonazzola, Phys. Rev. D 63, 064029 (2001). 
We set G = c = K = 1, where k is the polytropic constant 
of the EOS: P = Kpo- For scaling with k, see ref 0. 

M. D. Duez, P. Marronetti, S. L. Shapiro and 
T. W. Baumgarte, Phys. Rev. D 67, 024004 (2003). 
M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 
(1995). T. W. Baumgarte and S. L. Shapiro, Phys. Rev. 

D, 59, 024007 (1999). 

F. D. Swesty, E. Y. Wang and A. C. Calder, Astrophys. 
J. 541, 937 (2000). 

The time coordinate has been transformed into fractions 
of the initial orbital period. Note that the periods differ 
for binaries at different separation. 

J. C. Lombardi, F. A. Rasio and S. L. Shapiro, Phys. 
Rev. D 56, 3416 (1997). 

The total angular momentum like the gravitational mass 
is not strictly conserved due to the emission of gravita- 
tional waves. However, the loss rate of angular momen- 
tum for a typical orbit at these separations is of the order 
of 1% to 2%. Deviations exceeding these values therefore 
indicate numerical error. 

This figure presents the L2 norm of the constraint resid- 
uals. The details of these plots as well as the normaliza- 
tions used here can be found in |T3| . In this paper the mo- 
mentum constraint curve shows the average of the three 
spatial components. 
M. Miller, gr-qc/0305024] 

M. Miller and W. M. Suen, gr-qc/0301112 (2003), and 



5 



M. Miller, P. Gressman and W. M. Suen, gr-qc/0312030 present in the initial data, which leaves M and J nearly 

(2003). unchanged. 
[25] E.g., the loss of a small amount of spurious radiation 



